home *** CD-ROM | disk | FTP | other *** search
/ Libris Britannia 4 / science library(b).zip / science library(b) / DDJMAG / DDJ9212.ZIP / fortran.asc < prev    next >
Text File  |  1992-11-30  |  14KB  |  622 lines
  1. _WHY C++ WILL REPLACE FORTRANB_
  2. by Thomas Keffer
  3.  
  4.  
  5. [LISTING ONE]
  6.  
  7. // This class holds a reference count and a pointer to the vector data
  8. class DataBlock {
  9. private:
  10.   unsigned short   refs;          // Number of references
  11.   void*            array;         // Pointer to raw, untyped data
  12. public:
  13.   DataBlock(unsigned nelem, size_t elemsize)
  14.     {array = new char[nelem*elemsize]; refs = 1; }
  15.   ~DataBlock() { delete array; }
  16.   void             addReference()      {refs++;}
  17.   unsigned short   removeReference()   {return --refs;}
  18.   void*            data()              {return array;}
  19. };
  20.  
  21.  
  22. [LISTING TWO]
  23.  
  24. class DoubleVec {
  25.   DataBlock*       block;    // Pointer to the DataBlock
  26.   double*          begin;    // Start of data
  27.   unsigned         npts;     // Length of the vector
  28.   int              step;     // Stride length
  29.   DoubleVec(const DoubleVec&, int, unsigned, int); // For slices
  30. public:
  31.   DoubleVec(unsigned n, double val);   // n elements, initialized to val
  32.   DoubleVec(const DoubleVec& a);       // Copy constructor
  33.   ~DoubleVec();
  34.   DoubleVec        slice(int start, int stride, unsigned lgt) const;
  35.   unsigned         length() const           {return npts;}
  36.   int              stride() const           {return step;}
  37.   double&          operator()(int i) const  {return begin[i*step];}
  38.   DoubleVec&       operator=(const DoubleVec& v);
  39. };
  40. // Copy constructor:
  41. DoubleVec::DoubleVec(const DoubleVec& a)
  42. {
  43.   block = a.block;
  44.   block->addReference();
  45.   npts = a.npts;
  46.   begin = a.begin;
  47.   step = a.step;
  48. }
  49. // Destructor:
  50. DoubleVec::~DoubleVec()
  51. {
  52.   if (block->removeReference()==0)
  53.     delete block;
  54. }
  55.  
  56.  
  57. [LISTING THREE]
  58.  
  59. // Construct a vector that is a slice of an existing vector:
  60. DoubleVec::DoubleVec(const DoubleVec& v, int start, unsigned n, int str)
  61. {
  62.      block = v.block;
  63.      block->addReference();
  64.      npts = n;
  65.      begin = v.begin + start*v.stride;    // Add start points
  66.      stride = str*v.stride;               // Accumulate strides
  67. }
  68.  
  69.  
  70.  
  71.  
  72. [LISTING FOUR]
  73.  
  74. class DoubleMatrix : private DoubleVec {
  75.   unsigned ncols;  // Number of columns
  76.   unsigned nrows;  // Number of rows
  77. public:
  78.   DoubleMatrix();
  79.   DoubleMatrix(unsigned rows, unsigned cols, double initval);
  80.   DoubleMatrix(const DoubleMatrix& m);           // Reference to m will be made
  81.   unsigned         cols() const {return ncols;}
  82.   unsigned         rows() const {return nrows;}
  83.   DoubleVec        col(unsigned j) const           // Return col as slice
  84.     { return DoubleVec::slice(j*nrows, 1, nrows); }
  85.  
  86.   DoubleVec        row(unsigned i) const;          // Return row as slice
  87.     { return DoubleVec::slice(i, nrows, ncols); }
  88.   DoubleVec        diagonal() const;               // Return diagonal as slice
  89.   double&          operator()(int i, int j) const; // Subscripting
  90.  ...        // Other operators
  91. };
  92.  
  93.  
  94.  
  95.  
  96. [LISTING FIVE]
  97.  
  98. class LUDecomp : private DoubleMatrix {
  99.    int*                     permute;  // Row permutations
  100. public:
  101.   LUDecomp(const DoubleMatrix&);
  102.   ~LUDecomp() { delete permute; }
  103.   int                   isSingular() const;
  104.   friend double         determinant(const LUDecomp&);
  105.   friend DoubleMatrix   inverse(const LUDecomp&);
  106.   friend DoubleVec      solve(const LUDecomp&, const DoubleVec&);
  107. };
  108.  
  109.  
  110.  
  111. [LISTING SIX]
  112.  
  113. class  FFTServer {
  114.    unsigned         npts;
  115.    ComplexVec       theRoots;
  116.    void             calculateRoots();
  117. public:
  118.    FFTServer(unsigned n = 0) {setOrder(n);}
  119.    void             setOrder(unsigned n)
  120.      { npts = n; calculateRoots(); }
  121.    ComplexVec       fourier(const ComplexVec& v){
  122.      if(v.length() != npts) setOrder(v.length());
  123.      ComplexVec tran;
  124.      // ... (calculate transform, put it in tran)
  125.      return tran;
  126.    }
  127.  };
  128.  
  129.  
  130.  
  131. Example 1
  132.  
  133. ComplexVec b;           // Declare a complex vector
  134. cin >> b;               // Read it in
  135. FFTServer s;            // Allocate an FFT server
  136.                         // Calculate the transform of b:
  137. ComplexVec theTransform = s.fourier(b);
  138. cout << theTransform;   // Print it out
  139.  
  140.  
  141.  
  142. Example 2
  143.  
  144. (a)
  145.  
  146. DoubleVec a;         // Null vector - no elements at all, but can be resized
  147. DoubleVec b(8);      // 8 elements long, uninitialized
  148. DoubleVec c(8,1);    // 8 elements, initialized to 1.0
  149. DoubleVec d(8,1,2);  // 8 elements, initialized to 1, 3, 5, ...
  150.  
  151.  
  152. (b)
  153.  
  154. b = 2.0;      // Set all elements of b to 2
  155. b = c + d;    // Set b to the sum of c and d
  156. b *= 2;       // Multiply each element in b by 2.
  157.  
  158.  
  159. (c)
  160.  
  161. b[2] = 4.0;   // Set the 3'rd element of b to 4.0.
  162. c[1] = b[3];  // Set the 2'nd element of c to the 4'th element of b
  163.  
  164.  
  165. (d)
  166.  
  167.  
  168. DoubleVec a(10, 0);     // 10 element vector
  169. for(int i = 0; i<10; i+=2)
  170.    a[i] = 1;
  171.  
  172.  
  173.  
  174. Example 3
  175.  
  176. (a)
  177.  
  178. a.slice(0, 2) = 1;  // Starting with element 0; set every other element to 1
  179.  
  180.  
  181. (b)
  182.  
  183. class DoubleVec {
  184.  ...
  185. public:
  186.  ...
  187.  
  188.  operator         DoubleSlice();   // For type conversion
  189.  DoubleSlice      slice(int start, int step, unsigned N) const;
  190.  };
  191.  // The "helper class":
  192.  class DoubleSlice {
  193.      DoubleVec*       theVector;
  194.      int              startElement;
  195.      unsigned         sliceLength;
  196.      int              step;
  197.  public:
  198.      ...
  199.      friend DoubleVec operator+(const DoubleSlice&, const DoubleSlice&);
  200.      ...
  201.  };
  202.  
  203. (c)
  204.  
  205. DoubleVec b(10, 0), c(10, 1);
  206. DoubleVec d = b.slice(0, 2) + c.slice(1, 2);
  207.  
  208. (d)
  209.  
  210. DoubleVec g = b + c;    // DoubleVec to DoubleSlice type conversion occurs
  211.  
  212.  
  213. Example 4
  214.  
  215. (a)
  216.             
  217. DoubleVec     real(const ComplexVec&);
  218.  
  219. (b)            
  220.  
  221. ComplexVec a(10, 0);         // (0,0), (0,0), (0,0), ...
  222. real(a) = 1.0;               // (1,0), (1,0), (1,0), ...
  223.  
  224.  
  225. Example 5
  226.  
  227. (a)
  228.  
  229.  
  230. DoubleMatrix a(10, 10, 0);   // 10 by 10 matrix, initialized to zero
  231. a.row(3) = 1;                // Set row 3 to 1
  232. a.col(2) = a.col(4);         // Copy column 4 to column 2
  233.  
  234.  
  235.  
  236.  
  237. (b)
  238.  
  239.  
  240. DoubleMatrix I(10, 10, 0); // 10x10 initialized to 0
  241. I.diagonal() = 1;          // Create an identity matrix
  242.  
  243.  
  244. Example 6
  245. Example 6
  246.  
  247. (a)
  248.  
  249. DoubleMatrix a(10, 10);
  250. // ... (initialize a somehow)
  251. // Construct the LU decomposition of a:
  252. LUDecomp aLU(a);
  253.  
  254. // Now use it:
  255. double det = determinant(aLU);
  256. DoubleMatrix aInverse = inverse(aLU);
  257.  
  258. (b)
  259.  
  260. // 5 different sets of linear equations to be solved:
  261. DoubleVec b[5], x[5];
  262. // ... (set up the 5 vectors b and the 5 vectors x, each
  263. //      with 10 elements as per the matrix a above)
  264. for(int i = 0; i < 5; i++)
  265.    x[i] = solve(aLU, b[i]);
  266.  
  267. (c) [SET THIS IN HELV, NOT PRESTIGE]
  268.  
  269. a xi = bi;   i = 0, ..., 4.
  270.  
  271.  
  272. (d)
  273.  
  274. DoubleMatrix a(10, 10);
  275. // ... (initialize a)
  276. // Calculate the inverse directly from a.
  277. // A DoubleMatrix to LUDecomp type conversion takes place automatically:
  278. DoubleMatrix aInverse = inverse(a);
  279.  
  280.  
  281. (e)
  282.  
  283. DoubleMatrix aInverse = inverse(LUDecomp(a));
  284.  
  285.  
  286.  
  287. Example 7
  288.  
  289. ComplexVec timeVector(30);
  290. FFTServer aServer;      // Allocate a server
  291. // Will automagically reconfigure for a vector of length 30:
  292. ComplexVec spectrum = aServer.fourier(timeVector);
  293.  
  294.  
  295. Example 8
  296.  
  297. class Force;                                             // 1
  298. class String {                                           // 2
  299. public:
  300.   String(double length, double tension, double density); // 3
  301.   void setPoints(int nx);                                // 4
  302.   void timeStep(double dt, const Force& force);          // 5
  303.   DoubleVec displacement() const;                        // 6
  304. private:
  305.   DoubleVec u;                                           // 7
  306.   double    cSquared;
  307.   double    length;
  308. };
  309.  
  310.  
  311. Example 9
  312.  
  313. (a)
  314.  
  315.  
  316. class Force {
  317. public:
  318.    virtual DoubleVec   value() = 0;      // 1
  319. };
  320.  
  321.  
  322. (b)
  323.  
  324. class WindForceString : public Force {
  325. public:
  326.   WindForceString( String& string, double dragCoeff );
  327.   void                 setVelocity(double windspeed);
  328.   virtual DoubleVec    value();          // 1
  329. private:
  330.   String&               myString;        // The string we are tracking
  331.   double                wind;            // Present wind speed.
  332.   double                drag;
  333. };
  334.  
  335.  
  336.  
  337. (c)
  338.  
  339.             
  340. DoubleVec WindForceString::value()
  341. {
  342.     // Get present string displacement:
  343.     DoubleVec d = string.displacement();
  344.     return -drag*wind*wind*d;  // Some (bogus) calculation
  345. }
  346.  
  347.  
  348.  
  349.  
  350. Figure 1:
  351.  
  352.  
  353. C++                                    Fortran
  354.  
  355. #include <dvec.h>                             parameter (M=8000)
  356. #include <stdlib.h>                           double precision a(M), b(M)
  357. #include <iostream.h>                         double precision c(M)
  358.                                              
  359. main(int argc, char* argv[]) {                read(5,*) N, Niter
  360.                                              
  361.    unsigned N      = atoi(argv[1]);           write(6,1000)Niter, N
  362.    unsigned Niter  = atoi(argv[2]);    1000   format(i8,' iterations of ',
  363.                                              *      i8, ' points each.')
  364.    cout << Niter << "iterations of "              
  365.      << N << " points each.\n";               do 50 j=1,N
  366.                                               a(j) = 3
  367.     DoubleVec a(N, 3), b(N, 5);        50     b(j) = 5 
  368.                                               
  369.     while(Niter--)                            do 100 i=1,Niter 
  370.        DoubleVec c = a * b;                          do 100 j=1,N
  371. }                                      100           c(j) = a(j) * b(j)
  372.                                               stop
  373.                                               end
  374.  
  375.  
  376.  
  377. Figure 2
  378.  
  379.  
  380.  
  381. (a) 
  382.  
  383. #include <dgenfct.h>
  384. #include <rstream.h>
  385.  
  386.  
  387.  
  388.  
  389. class DTestMatrix : public DoubleGenMat {
  390. public:
  391.     DTestMatrix(unsigned order);
  392. };
  393.  
  394. class DTestRHS : public DoubleVec {
  395. public:
  396.     DTestRHS(const DTestMatrix&);
  397. };
  398.  
  399. double second();
  400. double epslon(double);
  401.  
  402.  
  403. const unsigned N = 90;
  404. const unsigned long ops = 2.0*N*N*N/3.0 + 2.0*N*N;
  405.  
  406. void main(){
  407.     DTestMatrix a(N);
  408.     double norma = maxVal(abs(a));
  409.     DTestRHS b(a);
  410.  
  411.     double t1 = second();
  412.     // Construct the LU Factorization:
  413.     DoubleGenFact fact(a);
  414.  
  415.     t1 = second() - t1;
  416.     double t2 = second();
  417.  
  418.     DoubleVec x = solve(fact, b);
  419.  
  420.     t2 = second() - t2;
  421.     double total = t1+t2;
  422.     double mflops = ops / (1.0e6*total);
  423.  
  424.     DoubleVec tol = a.product(x) - b;
  425.  
  426.  
  427.  
  428.  
  429.  
  430.  
  431.  
  432.  
  433.     double resid = maxVal(abs(tol));
  434.     double normx = maxVal(abs(x));
  435.  
  436.  
  437.  
  438.  
  439.  
  440.     double eps = epslon(1.0);
  441.     double residn= resid / (N*norma*normx*eps);
  442.  
  443.     cout << "Normalized residual    = " << residn << NL;
  444.     cout << "Residual               = " << resid  << NL;
  445.     cout << "Machine precision      = " << eps    << NL;
  446.     cout << "Factorization time     = " << t1     << NL;
  447.     cout << "Solution time          = " << t2     << NL;
  448.     cout << "Total time             = " << total  << NL;
  449.     cout << "MFLOPS                 = " << mflops << NL;
  450. }
  451.  
  452. DTestMatrix::DTestMatrix(unsigned n) : DoubleGenMat(n,n) {
  453.  
  454.  
  455.     long init = 1325;
  456.     for(int j=0; j<n; j++){
  457.       for(int i=0; i<n; i++){
  458.         init = 3125*init % 65536;
  459.         sub(i,j) = (init-32768.0)/16384.0;
  460.       }
  461.     }
  462. }
  463.  
  464.  
  465. DTestRHS::DTestRHS(const DTestMatrix& a) : 
  466.     DoubleVec(a.rows(), 0.0) {
  467.  
  468.  
  469.  
  470.     for(int i=0; i<length(); i++)
  471.       (*this)(i) = sum(a.row(i));
  472. }
  473.  
  474.  
  475.  
  476.  
  477.  
  478. Total (nonblank) lines  = 52
  479. Executable size = 64868 bytes
  480.   (Borland C++ V2.0 w. optimizer & 8087, large memory model)
  481.  
  482. 16 MHz 386 w. 80387:
  483.  
  484. Normalized residual = 1.24482
  485. Residual            =  .4973799e-13
  486. Machine precision   =  .2220446e-15
  487. Factorization time  = 3.97
  488. Solution time       = 0.13
  489. Total time          = 4.10
  490. MFLOPS              = 0.122
  491.  
  492.  
  493.  
  494.  
  495.  
  496.  
  497. (b) 
  498.  
  499.      double precision a(90,90),b(90),x(90)
  500.      double precision ops, mflops, norma, normx
  501.      double precision resid, residn, eps
  502.      double precision t1, t2, total
  503.      integer ipvt(90)
  504.  
  505.  
  506.  
  507.  
  508.  
  509.  
  510.  
  511.  
  512.  
  513.  
  514.  
  515.      double precision second
  516.      double precision epslon
  517.  
  518.      lda = 90
  519.      n = 90
  520.      ops = (2.0e0*n**3)/3.0e0 + 2.0e0*n**2
  521.  
  522.  
  523.      call matgen(a,lda,n,b,norma)
  524.  
  525.  
  526.  
  527.      t1 = second()
  528.  
  529.      call dgefa(a,lda,n,ipvt,info)
  530.  
  531.      t1 = second() -t1
  532.      t2 = second()
  533.  
  534.      call dgesl(a,lda,n,ipvt,b,0)
  535.  
  536.      t2 = second() - t2
  537.      total = t1+t2
  538.      mflops = ops/(1.0e6*total)
  539.  
  540.  
  541.  
  542.  
  543.      do 10 i = 1,n
  544.        x(i) = b(i)
  545.   10 continue
  546.      call matgen(a,lda,n,b,norma)
  547.      do 20 i = 1,n
  548.        b(i) = -b(i)
  549.   20 continue
  550.      CALL DMXPY(n,b,n,lda,x,a)
  551.  
  552.      RESID = 0.0
  553.      NORMX = 0.0
  554.      DO 30 I = 1,N
  555.        RESID = amax1( RESID, ABS(b(i)) )
  556.        NORMX = amax1( NORMX, ABS(X(I)) )
  557.   30 CONTINUE
  558.  
  559.      eps = epslon(1.0D0)
  560.      RESIDn = RESID/( N*NORMA*NORMX*EPS )
  561.  
  562.      write(6,1000)RESIDn
  563. 1000 format(' Normalized residual    =', g16.7)
  564.      write(6,1001)RESID
  565. 100  format(' Residual               = ', g16.7)
  566.      write(6,1002)eps
  567. 1002 format(' Machine precision      = ', g16.7)
  568.      write(6,1003)t1
  569. 1003 format(' Factorization time     = ', g16.7)
  570.      write(6,1004)t2
  571. 1004 format(' Solution time          = ', g16.7)
  572.      write(6,1005)total
  573. 1005 format(' Total time             = ', g16.7)
  574.      write(6,1006)mflops
  575. 1006 format(' MFLOPS                 = ', g16.7)
  576.      stop
  577.      end
  578.  
  579.  
  580.      subroutine matgen(a,lda,n,b,norma)
  581.      double precision a(lda,1),b(1),norma
  582.  
  583.      init = 1325
  584.      norma = 0.0
  585.      do 30 j = 1,n
  586.        do 20 i = 1,n
  587.          init = mod(3125*init,65536)
  588.          a(i,j) = (init - 32768.0)/16384.0
  589.          norma = amax1(a(i,j), norma)
  590.   20   continue
  591.   30 continue
  592.  
  593.      do 35 i = 1,n
  594.        b(i) = 0.0
  595.   35 continue
  596.  
  597.      do 50 j = 1,n
  598.        do 40 i = 1,n
  599.          b(i) = b(i) + a(i,j)
  600.   40   continue
  601.   50 continue
  602.      return
  603.      end
  604.  
  605.  
  606. Total (nonblank) lines  = 71
  607. Executable size = 87482 bytes
  608.   (Microsoft Fortran V5.0 w. optimizer & 8087)
  609.  
  610. 16 MHz 386 w. 80387:
  611.  
  612. Normalized residual  =     1.212636    
  613. Residual             =      .4843265E-13
  614. Machine precision    =      .2220446E-15
  615. Factorization time   =     4.230000    
  616. Solution time        =      .170000    
  617. Total time           =     4.400000    
  618. MFLOPS               =      .1141364    
  619.  
  620.        
  621.  
  622.